Geo-preprocessing

Set Up

# Packages ---------------------------------------------------------------
# (uncomment next 2 lines if you want auto-install)
# pkgs <- c("leaflet","terra","sf","leafem","raster")
# to_install <- pkgs[!sapply(pkgs, requireNamespace, quietly = TRUE)]; if (length(to_install)) install.packages(to_install)

library(terra)    # raster
library(sf)       # vectors
library(leaflet)  # web map
library(leafem)   # helpers for raster in leaflet
library(raster)   # addRasterImage prefers RasterLayer
library(igraph)
library(units)
library(ggplot2)
library(smoothr)
library(dplyr)

Read Data

coralSDM <- terra::rast("data/maxent_predrast_GBR_AhyaD_015_lq2.tif")
crs(coralSDM)

[1] "PROJCRS[\"GDA2020 / MGA zone 55\",\n    BASEGEOGCRS[\"GDA2020\",\n        DATUM[\"Geocentric Datum of Australia 2020\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",7844]],\n    CONVERSION[\"Map Grid of Australia zone 55\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",147,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",10000000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"Australia - onshore and offshore between 144°E and 150°E.\"],\n        BBOX[-50.89,144,-9.23,150.01]],\n    ID[\"EPSG\",7855]]"

bb_wgs <- ext(147.0, 147.1, -18.66, -18.6)

# 2) Turn extent -> polygon with CRS, then project to raster CRS (MGA55)
bb_wgs_poly <- as.polygons(bb_wgs, crs = "EPSG:4326")
bb_mga_poly <- project(bb_wgs_poly, crs(coralSDM))

# 3) Crop (and optionally mask) using on-disk processing (good for big rasters)
coralSDM_crop <- crop(coralSDM, bb_mga_poly)

# Optional (keeps only inside the polygon edges):
# out <- mask(out, bb_mga_poly, filename = "outputs/coralSDM_clip.tif", overwrite = TRUE)


plot(coralSDM_crop)

Determine Tow Path

#-----------------------------
# 1) Load raster (probability 0-1)
#-----------------------------
r <- coralSDM_crop  # CRS = GDA2020 / MGA 55 (metres)
crs_r <- st_crs(crs(r))  # parse the WKT from terra into an sf crs


#-----------------------------
# 2) Find probability thresholds by *cumulative probability mass*
#   t50: minimal value so that pixels >= t50 sum to 50% of total mass
#   t85: pixels >= t85 sum to 85% of total mass (gives "next ~35%" between t85 and t50)
#-----------------------------
v <- values(r, mat = FALSE, na.rm = TRUE)
o <- sort(v, decreasing = TRUE)
cs <- cumsum(o) / sum(o)

t50 <- o[which(cs >= 0.50)[1]]
t85 <- o[which(cs >= 0.85)[1]]  # change 0.85 to 0.80 or 0.90 if you want 30% or 40%

# Masks
top50     <- r >= t50
next35ish <- (r >= t85) & (r < t50)

# make them 1/NA so polygons only represent the TRUE area
top50     <- ifel(top50, 1, NA)
next35ish <- ifel(next35ish, 1, NA)
#-----------------------------
# 3) Raster->polygon (dissolved) for each zone
#-----------------------------
# polygonise; no need to filter on a field name
poly_top50     <- as.polygons(top50,     dissolve = TRUE)      |> st_as_sf()
poly_next35ish <- as.polygons(next35ish, dissolve = TRUE)      |> st_as_sf()




# Convert to sf for path generation
sf_top50     <- st_as_sf(poly_top50)
sf_next35ish <- st_as_sf(poly_next35ish)

# If sf_top50 / sf_next35ish came from terra::as.polygons(), still force-transform:
sf_top50     <- st_transform(sf_top50,     crs_r)
sf_next35ish <- st_transform(sf_next35ish, crs_r)

# Optional: clean tiny slivers
drop_small <- function(sf_obj, min_area_m2 = 2000) {
  sf_obj %>% mutate(a = st_area(geometry)) %>% filter(a >= set_units(min_area_m2, m^2)) %>% st_make_valid()
}
sf_top50     <- drop_small(sf_top50,     2000)
sf_next35ish <- drop_small(sf_next35ish, 2000)

#-----------------------------
# 4) Make a boustrophedon tow plan inside a polygon
#   spacing_m: distance between adjacent lines (≈ swath width)
#   angle_deg: optional rotation to align with reef orientation
#-----------------------------


make_tow_path <- function(poly_sf, spacing_m = 75, angle_deg = 0) {
  stopifnot(inherits(poly_sf, "sf"))
  if (isTRUE(st_is_longlat(poly_sf)[1])) stop("Polygon must be in a projected CRS (metres).")

  poly_sf <- st_make_valid(poly_sf)
  crs_    <- st_crs(poly_sf)

  # Single unified polygon (sfc)
  poly_u <- st_union(poly_sf)
  ctr    <- st_centroid(poly_u)

  # Rotation matrices
  rad  <- angle_deg * pi/180
  R    <- matrix(c(cos(rad), -sin(rad), sin(rad), cos(rad)),  nrow = 2)
  Rinv <- matrix(c(cos(-rad), -sin(-rad), sin(-rad), cos(-rad)), nrow = 2)

  # Rotate polygon; IMPORTANT: re-attach CRS after affine math
  g <- (poly_u - ctr) * R + ctr
  st_crs(g) <- crs_
  poly_r <- st_sf(geometry = g)

  # Build vertical lines across (slightly extended) rotated bbox
  bb   <- st_bbox(poly_r)
  span <- max(bb["xmax"] - bb["xmin"], bb["ymax"] - bb["ymin"])
  margin <- 0.1 * span
  xs <- seq(bb["xmin"] - margin, bb["xmax"] + margin, by = spacing_m)

  segs <- lapply(xs, function(x)
    st_linestring(matrix(c(x, bb["ymin"] - margin, x, bb["ymax"] + margin),
                         ncol = 2, byrow = TRUE))
  )
  lines <- st_sf(geometry = st_sfc(segs, crs = crs_))

  # Intersect (same CRS on both)
  stopifnot(identical(st_crs(lines), st_crs(poly_r)))
  tran <- suppressWarnings(st_intersection(lines, poly_r))
  if (!nrow(tran)) stop("No transects created. Check spacing/angle/extent.")

  # Order left->right and snake
  tran$xx <- st_coordinates(st_centroid(tran))[, 1]
  tran <- tran[order(tran$xx), ]

  geoms <- lapply(seq_len(nrow(tran)), function(i){
    coords <- st_coordinates(tran[i, ])[, 1:2, drop = FALSE]
    if (nrow(coords) < 2) return(st_linestring(coords))
    if (i %% 2 == 0) coords <- coords[nrow(coords):1, , drop = FALSE]
    st_linestring(coords)
  })
  st_geometry(tran) <- st_sfc(geoms, crs = crs_)

  # Merge segments, rotate back, re-attach CRS
  path_rot <- st_line_merge(st_union(tran))
  if (inherits(path_rot, "sfc_GEOMETRYCOLLECTION"))
    path_rot <- st_collection_extract(path_rot, "LINESTRING")
  st_crs(path_rot) <- crs_

  g2 <- (path_rot - ctr) * Rinv + ctr
  st_crs(g2) <- crs_
  st_sf(geometry = g2)
}




# Build paths (tune spacing/angle as needed)
tow_top50     <- make_tow_path(sf_top50,     spacing_m = 75, angle_deg = 20)  # "fast" survey
tow_next35ish <- make_tow_path(sf_next35ish, spacing_m = 75, angle_deg = 20)  # "rigorous" extension

#-----------------------------
# 5) Save results
#-----------------------------
st_write(sf_top50,     "outputs/top50_poly.gpkg",     delete_dsn = TRUE)
st_write(sf_next35ish, "outputs/next35_poly.gpkg",    delete_dsn = TRUE)
st_write(tow_top50,     "outputs/tow_top50.gpkg",     delete_dsn = TRUE)
st_write(tow_next35ish, "outputs/tow_next35.gpkg",    delete_dsn = TRUE)

# --- INPUTS ---
r <- coralSDM_crop   # SpatRaster of probability 0..1 in MGA55 (metres)
summary(values(r))

     lyr.1          
 Min.   :0.0002913  
 1st Qu.:0.0053625  
 Median :0.0103037  
 Mean   :0.0355549  
 3rd Qu.:0.0275046  
 Max.   :0.9999997  
 NA's   :16         

# --- Build CONDUCTANCE from probability (safer than cost=1-p) ---
alpha=1.5
cond  <- app(r, function(x) (pmax(0, pmin(1, x)))^alpha)
cond[is.na(cond)] <- 0

plot(cond)

hist(values(cond), main = "Conductance Values", xlab = "Conductance")

# gdistance needs RasterLayer
cond_r <- raster::raster(cond)

# Transition (use mean because we're using conductance)
tr <- gdistance::transition(cond_r, function(x) mean(x, na.rm = TRUE), directions = 8)
tr <- gdistance::geoCorrection(tr, type = "c")

# Get cells with high conductance (you can lower threshold if needed)
high_cells <- which(raster::values(cond_r) > 0.8)
length(high_cells)  # You want at least 2!

[1] 1040

set.seed(42)  # for reproducibility
chosen_cells <- sample(high_cells, 2)
chosen_coords <- raster::xyFromCell(cond_r, chosen_cells)

# Assign as new start/end coordinates
start_xy <- chosen_coords[1, ]
end_xy   <- chosen_coords[2, ]


# Start/end as SpatialPoints (CRS preserved)
crs_ <- st_crs(crs(r))                       # parse terra WKT to sf crs
p1   <- st_sfc(st_point(start_xy), crs = crs_) |> as("Spatial")
p2   <- st_sfc(st_point(end_xy),   crs = crs_) |> as("Spatial")

snap_to_tr <- function(pt_sp, tr, rRL) {
  stopifnot(inherits(pt_sp, "SpatialPoints"))

  # Use raster to get valid cell indices (non-zero conductance)
  valid_cells <- which(raster::values(rRL) > 0)

  if (!length(valid_cells)) stop("No valid conductance cells in raster.")

  # Check if point is already valid
  cell0 <- raster::cellFromXY(rRL, sp::coordinates(pt_sp))
  ok <- is.finite(cell0) && (cell0 %in% valid_cells)

  if (!ok) {
    # Snap to nearest valid cell center
    valid_xy <- raster::xyFromCell(rRL, valid_cells)
    xy <- sp::coordinates(pt_sp)
    i  <- which.min((valid_xy[,1] - xy[1,1])^2 + (valid_xy[,2] - xy[1,2])^2)
    pt_sp <- SpatialPoints(matrix(valid_xy[i,], ncol = 2), proj4string = raster::crs(rRL))
  }

  pt_sp
}
# 
# 
# # --- make sure p1/p2 are in the correct CRS and snapped to the graph ---
# # If your points started as sf/lonlat, reproject first:
# # pts_sf <- st_transform(pts_sf_ll, st_crs(crs(r)))
# # p1 <- as(pts_sf[1,], "Spatial"); p2 <- as(pts_sf[2,], "Spatial")
# 
# p1 <- snap_to_tr(p1, tr, cond_r)
# p2 <- snap_to_tr(p2, tr, cond_r)

# Shortest path (as SpatialLines), then to sf
sPath <- gdistance::shortestPath(tr, p1, p2, output = "SpatialLines")
sPath_sf <- sf::st_as_sf(sPath)

# --- Plot with ggplot (convert raster to df once) ---
df <- as.data.frame(r, xy = TRUE, na.rm = FALSE)
names(df)[3] <- "p"

ggplot(df) +
  geom_raster(aes(x = x, y = y, fill = p)) +
  geom_sf(data = sPath_sf, colour = "orange", linewidth = 2) +
  # coord_sf(crs = crs_$wkt) +
  scale_fill_viridis_c(name = "p(coral)", limits = c(0,1), oob = scales::squish) +
  theme_bw()

# Convert raster to binary mask (1 = high coral prob)
reef_mask <- cond > 0.8

# Convert to polygons (this will outline the reef)
reef_poly <- as.polygons(reef_mask, dissolve = TRUE) |> st_as_sf() |>
  dplyr::filter(lyr.1==1)

# Get raster extent as an sf polygon
r_extent <- st_as_sfc(st_bbox(cond))

# Simplify the largest ring

reef_ring <- reef_poly[which.max(st_area(reef_poly)), ]  # take the biggest one
# Id otn think above is doing anythin usefuk as the polygons are just a single mulitplygon

reef_line <- reef_poly |>
  st_cast("MULTILINESTRING") |>
  st_cast("LINESTRING")

# 2. Sample points along the boundary (e.g., 20)
pts <- st_line_sample(reef_line, n = 1, type = "regular") |> 
  st_cast("POINT") |> as("Spatial")


# 4. Snap all points to graph
# pts_snapped <- lapply(pts_sp, snap_to_tr, tr = tr, rRL = cond_r)
pts_snapped <- lapply(seq_len(length(pts)), function(i) {
  snap_to_tr(pts[i, ], tr, cond_r)
})

# 5. Loop through pairs to build full path
paths_list <- list()
for (i in seq_along(pts_snapped)[-length(pts_snapped)]) {
  print(i)
  paths_list[[i]] <- gdistance::shortestPath(tr, pts_snapped[[i]], pts_snapped[[i+1]], output = "SpatialLines")
}

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12
[1] 13
[1] 14
[1] 15
[1] 16
[1] 17
[1] 18
[1] 19
[1] 20
[1] 21
[1] 22
[1] 23
[1] 24
[1] 25
[1] 26
[1] 27
[1] 28
[1] 29
[1] 30
[1] 31
[1] 32
[1] 33
[1] 34
[1] 35
[1] 36
[1] 37
[1] 38
[1] 39
[1] 40
[1] 41
[1] 42
[1] 43
[1] 44
[1] 45
[1] 46
[1] 47
[1] 48
[1] 49
[1] 50
[1] 51
[1] 52
[1] 53
[1] 54
[1] 55
[1] 56
[1] 57
[1] 58
[1] 59
[1] 60
[1] 61
[1] 62
[1] 63
[1] 64
[1] 65
[1] 66
[1] 67
[1] 68
[1] 69
[1] 70
[1] 71
[1] 72
[1] 73
[1] 74
[1] 75
[1] 76
[1] 77
[1] 78
[1] 79
[1] 80
[1] 81
[1] 82
[1] 83
[1] 84
[1] 85
[1] 86
[1] 87
[1] 88
[1] 89

# Optional: Close the loop
paths_list[[length(pts_snapped)]] <- gdistance::shortestPath(tr, pts_snapped[[length(pts_snapped)]], pts_snapped[[1]], output = "SpatialLines")

# 6. Combine into one sf object
paths_sf <- do.call(rbind, lapply(paths_list, sf::st_as_sf))

# Get CRS from your raster
reef_crs <- crs(coralSDM_crop, proj = TRUE)  # terra crs as WKT

# Assign to paths_sf
st_crs(paths_sf) <- reef_crs

df <- as.data.frame(r, xy = TRUE, na.rm = FALSE)
names(df)[3] <- "p"

ggplot(df) +
  geom_raster(aes(x = x, y = y, fill = p)) +
  geom_sf(data = paths_sf, colour = "orange", linewidth = 1) +
  # coord_sf(crs = crs_$wkt) +
  scale_fill_viridis_c(name = "p(coral)", limits = c(0,1), oob = scales::squish) +
  theme_bw()

Leaflet Map

rng <- as.numeric(terra::global(coralSDM_crop, "range", na.rm = TRUE))
pal_coral <- colorNumeric("magma", rng, na.color = "transparent")
paths_sf <- st_transform(paths_sf, 4326)
m <- leaflet() |>
  addProviderTiles(providers$Esri.OceanBasemap, group = "Esri Ocean") |>
  addProviderTiles(providers$Esri.WorldImagery, group = "Esri Satellite") |>
  addProviderTiles(providers$CartoDB.Positron,  group = "Light") |>
  addRasterImage(
    coralSDM_crop,
    colors  = pal_coral,
    opacity = 0.7,
    project = TRUE,
    group   = "CoralSDM"
  ) |>
  addRasterImage(
    cond_r,
    opacity = 0.7,
    project = TRUE,
    group   = "Conductance"
  ) |>
    addPolylines(
    data  = paths_sf,
    color = "springgreen",
    weight = 3,
    group = "Tow Path"
  ) |>
  addLegend(
    pal = pal_coral, values = rng, title = "Probability",
    position = "bottomright", opacity = 1
  ) |>

  addLayersControl(
  baseGroups    = c("Esri Ocean", "Esri Satellite", "Light"),
  overlayGroups = c("CoralSDM", "Conductance", "Tow Path"),
  options = layersControlOptions(collapsed = FALSE)
)|>
  leafem::addMouseCoordinates()

m

50th Percentile Paths

#-----------------------------
# 1) Load raster (probability 0-1)
#-----------------------------
r <- coralSDM_crop  # CRS = GDA2020 / MGA 55 (metres)
crs_r <- st_crs(crs(r))  # parse the WKT from terra into an sf crs

# Downsample to coarser resolution (e.g., factor of 2 or 5)
r <- terra::aggregate(coralSDM_crop, fact = 3, fun = mean)

#-----------------------------
# 2) Find probability thresholds by *cumulative probability mass*
#   t50: minimal value so that pixels >= t50 sum to 50% of total mass
#   t85: pixels >= t85 sum to 85% of total mass (gives "next ~35%" between t85 and t50)
#-----------------------------
v <- values(r, mat = FALSE, na.rm = TRUE)
v <- ifelse(v>0.2, v, NA)
o <- sort(v, decreasing = TRUE)
cs <- cumsum(o) / sum(o)

t50 <- o[which(cs >= 0.50)[1]]
 # change 0.85 to 0.80 or 0.90 if you want 30% or 40%

# Masks
top50     <- r >= t50


# make them 1/NA so polygons only represent the TRUE area
top50     <- ifel(top50, 1, NA)

#-----------------------------
# 3) Raster->polygon (dissolved) for each zone
#-----------------------------
# polygonise; no need to filter on a field name
poly_top50     <- as.polygons(top50,     dissolve = TRUE)      |> st_as_sf()


# Convert to sf for path generation
sf_top50     <- st_as_sf(poly_top50)

# If sf_top50 / sf_next35ish came from terra::as.polygons(), still force-transform:
sf_top50     <- st_transform(sf_top50,     crs_r)

# Optional: clean tiny slivers
drop_small <- function(sf_obj, min_area_m2 = 2000) {
  sf_obj %>% mutate(a = st_area(geometry)) %>% filter(a >= set_units(min_area_m2, m^2)) %>% st_make_valid()
}
sf_top50     <- drop_small(sf_top50,     2000)

reef_line <- sf_top50 |>
  # st_simplify(dTolerance = 100) |>
  st_cast("MULTILINESTRING") |>
  st_cast("LINESTRING")

# 2. Sample points along the boundary (e.g., 20)
pts <- st_line_sample(reef_line, n = 1, type = "regular") |> 
  st_cast("POINT") |> as("Spatial")


# 4. Snap all points to graph
# pts_snapped <- lapply(pts_sp, snap_to_tr, tr = tr, rRL = cond_r)
pts_snapped <- lapply(seq_len(length(pts)), function(i) {
  snap_to_tr(pts[i, ], tr, cond_r)
})

# Convert snapped SpatialPoints to sf
pts_sf <- st_as_sf(do.call(rbind, pts_snapped))

# Ensure CRS is set
st_crs(pts_sf) <- st_crs(reef_crs)

# Compute centroid of the entire reef (or polygon)
reef_centroid <- st_centroid(st_union(reef_line))

# Extract coordinates
coords <- st_coordinates(pts_sf)
centroid_coords <- st_coordinates(reef_centroid)

# Calculate angle from centroid to each point
angles <- atan2(coords[, 2] - centroid_coords[2], coords[, 1] - centroid_coords[1])

# Add angle as column
pts_sf$angle <- angles

# Sort clockwise
pts_sf_sorted <- pts_sf[order(angles), ]

# Replace snapped list with sorted order
pts_snapped <- lapply(1:nrow(pts_sf_sorted), function(i) {
  as(st_geometry(pts_sf_sorted[i, ]), "Spatial")
})


# 5. Loop through pairs to build full path
paths_list <- list()
for (i in seq_along(pts_snapped)[-length(pts_snapped)]) {
  print(i)
  paths_list[[i]] <- gdistance::shortestPath(tr, pts_snapped[[i]], pts_snapped[[i+1]], output = "SpatialLines")
}

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12
[1] 13
[1] 14
[1] 15
[1] 16
[1] 17
[1] 18
[1] 19
[1] 20
[1] 21
[1] 22
[1] 23
[1] 24
[1] 25
[1] 26
[1] 27
[1] 28
[1] 29
[1] 30
[1] 31
[1] 32
[1] 33
[1] 34
[1] 35
[1] 36
[1] 37
[1] 38
[1] 39
[1] 40
[1] 41
[1] 42
[1] 43

# Optional: Close the loop
paths_list[[length(pts_snapped)]] <- gdistance::shortestPath(tr, pts_snapped[[length(pts_snapped)]], pts_snapped[[1]], output = "SpatialLines")

# 6. Combine into one sf object
paths_sf <- do.call(rbind, lapply(paths_list, sf::st_as_sf))

# Get CRS from your raster
reef_crs <- crs(coralSDM_crop, proj = TRUE)  # terra crs as WKT

# Assign to paths_sf
st_crs(paths_sf) <- reef_crs

df <- as.data.frame(r, xy = TRUE, na.rm = FALSE)
names(df)[3] <- "p"

# 1. Extract coordinates from each SpatialLines path
coords_list <- lapply(paths_list, function(x) {
  if (!is.null(x)) {
    sp::coordinates(x)[[1]][[1]]  # extract path coords as matrix
  } else {
    NULL
  }
})
# 2. Filter out NULLs (e.g., failed paths)
coords_list <- Filter(Negate(is.null), coords_list)

# Combine all coordinate matrices into one continuous path
all_coords <- do.call(rbind, coords_list)

# Create a LINESTRING
path_line <- st_linestring(all_coords)

# Wrap into sf object
path_line_sf <- st_sf(geometry = st_sfc(path_line), crs = reef_crs)


path_smooth <- smooth(path_line_sf, method = "ksmooth", smoothness = 10)

ggplot(df) +
  geom_raster(aes(x = x, y = y, fill = p)) +
  geom_sf(data = paths_sf, colour = "orange", linewidth = 1) +
  geom_sf(data = path_smooth, colour = "lightgreen", linewidth = 1) +
  # coord_sf(crs = crs_$wkt) +
  scale_fill_viridis_c(name = "p(coral)", limits = c(0,1), oob = scales::squish) +
  theme_bw()

Create buffer for next 30% of habitat

# Downsample to coarser resolution (e.g., factor of 2 or 5)
r <- terra::aggregate(coralSDM_crop, fact = 3, fun = mean)
# 1. Buffer the first path (in meters, assuming projected CRS)
path_50 <- path_line_sf
path_50 <- st_transform(path_50, crs(r))

path_50_buffer <- sf::st_buffer(path_50, dist = 200)

# 2. Rasterize the buffer (must match resolution & extent of conductance raster)
buffer_raster <- terra::rasterize(path_50_buffer, r)
buffer_raster_aligned <- terra::project(buffer_raster, r, method = "near")
# 3. Mask the conductance raster to exclude the buffered path
r.masked <- terra::mask(r, buffer_raster_aligned, maskvalues = 1, inverse = F)
plot(r.masked)

buffer_raster_aligned <- terra::project(buffer_raster, cond, method = "near")
cond_30_masked <- terra::mask(cond, buffer_raster_aligned, maskvalues = 1, inverse = F)
plot(cond_30_masked)

# 3. Convert to RasterLayer (gdistance doesn't work with SpatRaster)
cond_30_r <- raster::raster(cond_30_masked)

# 4. Build a NEW transition matrix (based on updated conductance)
tr2 <- gdistance::transition(cond_30_r, function(x) mean(x, na.rm = TRUE), directions = 8)
tr2 <- gdistance::geoCorrection(tr2, type = "c")
### next 35

# 1) Load raster (probability 0-1)
#-----------------------------
r <- r.masked  # CRS = GDA2020 / MGA 55 (metres)
crs_r <- st_crs(crs(r))  # parse the WKT from terra into an sf crs

#-----------------------------
# 2) Find probability thresholds by *cumulative probability mass*
#   t50: minimal value so that pixels >= t50 sum to 50% of total mass
#   t85: pixels >= t85 sum to 85% of total mass (gives "next ~35%" between t85 and t50)
#-----------------------------
v <- values(r, mat = FALSE, na.rm = TRUE)
v <- ifelse(v>0.1, v, NA)
o <- sort(v, decreasing = TRUE)
cs <- cumsum(o) / sum(o)

t50 <- o[which(cs >= 0.50)[1]]
 # change 0.85 to 0.80 or 0.90 if you want 30% or 40%

# Masks
top50 <- r >= t50

# make them 1/NA so polygons only represent the TRUE area
top50 <- ifel(top50, 1, NA)

#-----------------------------
# 3) Raster->polygon (dissolved) for each zone
#-----------------------------
# polygonise; no need to filter on a field name
poly_top50 <- as.polygons(top50, dissolve = TRUE) |> st_as_sf()


# Convert to sf for path generation
sf_top50 <- st_as_sf(poly_top50)

# If sf_top50 / sf_next35ish came from terra::as.polygons(), still force-transform:
sf_top50 <- st_transform(sf_top50,crs_r)

# Optional: clean tiny slivers
drop_small <- function(sf_obj, min_area_m2 = 2000) {
  sf_obj %>% mutate(a = st_area(geometry)) %>% filter(a >= set_units(min_area_m2, m^2)) %>% st_make_valid()
}
sf_top50     <- drop_small(sf_top50,     2000)

reef_line.next <- sf_top50 |>
  # st_simplify(dTolerance = 100) |>
  st_cast("MULTILINESTRING") |>
  st_cast("LINESTRING")

# 2. Sample points along the boundary (e.g., 20)
pts.next <- st_line_sample(reef_line.next, n = 1, type = "regular") |> 
  st_cast("POINT") |> as("Spatial")


# 4. Snap all points to graph
# pts_snapped <- lapply(pts_sp, snap_to_tr, tr = tr, rRL = cond_r)
pts_snapped.next <- lapply(seq_len(length(pts.next)), function(i) {
  snap_to_tr(pts.next[i, ], tr2, cond_30_r)
})

# Convert snapped SpatialPoints to sf
pts_sf.next <- st_as_sf(do.call(rbind, pts_snapped.next))

# Ensure CRS is set
st_crs(pts_sf.next) <- st_crs(reef_crs)

# Compute centroid of the entire reef (or polygon)
reef_centroid <- st_centroid(st_union(reef_line))

# Extract coordinates
coords <- st_coordinates(pts_sf.next)
centroid_coords <- st_coordinates(reef_centroid)

# Calculate angle from centroid to each point
angles <- atan2(coords[, 2] - centroid_coords[2], coords[, 1] - centroid_coords[1])

# Add angle as column
pts_sf.next$angle <- angles

# Sort clockwise
pts_sf_sorted.next <- pts_sf.next[order(angles), ]

# Replace snapped list with sorted order
pts_snapped.next <- lapply(1:nrow(pts_sf_sorted.next), function(i) {
  as(st_geometry(pts_sf_sorted.next[i, ]), "Spatial")
})


# 5. Loop through pairs to build full path
paths_list.next <- list()
for (i in seq_along(pts_snapped.next)[-length(pts_snapped.next)]) {
  print(i)
  paths_list.next[[i]] <- gdistance::shortestPath(tr, pts_snapped.next[[i]], pts_snapped.next[[i+1]], output = "SpatialLines")
}

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12
[1] 13
[1] 14
[1] 15
[1] 16
[1] 17
[1] 18
[1] 19
[1] 20
[1] 21
[1] 22
[1] 23
[1] 24
[1] 25
[1] 26
[1] 27
[1] 28
[1] 29
[1] 30
[1] 31
[1] 32
[1] 33
[1] 34
[1] 35
[1] 36
[1] 37
[1] 38
[1] 39
[1] 40
[1] 41
[1] 42
[1] 43
[1] 44
[1] 45
[1] 46
[1] 47
[1] 48
[1] 49
[1] 50
[1] 51
[1] 52
[1] 53
[1] 54
[1] 55
[1] 56
[1] 57
[1] 58
[1] 59
[1] 60
[1] 61
[1] 62
[1] 63
[1] 64
[1] 65
[1] 66
[1] 67
[1] 68
[1] 69
[1] 70
[1] 71
[1] 72
[1] 73
[1] 74
[1] 75
[1] 76
[1] 77
[1] 78
[1] 79
[1] 80
[1] 81
[1] 82
[1] 83
[1] 84
[1] 85
[1] 86
[1] 87
[1] 88
[1] 89
[1] 90

# Optional: Close the loop
paths_list.next[[length(pts_snapped.next)]] <- gdistance::shortestPath(tr, pts_snapped.next[[length(pts_snapped.next)]], pts_snapped.next[[1]], output = "SpatialLines")

# 6. Combine into one sf object
paths_sf.next <- do.call(rbind, lapply(paths_list.next, sf::st_as_sf))

# Assign to paths_sf
st_crs(paths_sf.next) <- reef_crs

# 1. Extract coordinates from each SpatialLines path
coords_list.next <- lapply(paths_list.next, function(x) {
  if (!is.null(x)) {
    sp::coordinates(x)[[1]][[1]]  # extract path coords as matrix
  } else {
    NULL
  }
})
# 2. Filter out NULLs (e.g., failed paths)
coords_list.next <- Filter(Negate(is.null), coords_list.next)

# Combine all coordinate matrices into one continuous path
all_coords.next <- do.call(rbind, coords_list.next)

# Create a LINESTRING
path_line.next <- st_linestring(all_coords.next)

# Wrap into sf object
path_line_sf.next <- st_sf(geometry = st_sfc(path_line.next), crs = reef_crs)


path_smooth.next <- smooth(path_line_sf.next, method = "ksmooth", smoothness = 10)

ggplot(df) +
  geom_raster(aes(x = x, y = y, fill = p)) +
  geom_sf(data = paths_sf, colour = "orange", linewidth = 1) +
  geom_sf(data = path_smooth, colour = "lightgreen", linewidth = 1) +
  geom_sf(data = path_smooth.next, colour = "purple", linewidth = 1) +
  # coord_sf(crs = crs_$wkt) +
  scale_fill_viridis_c(name = "p(coral)", limits = c(0,1), oob = scales::squish) +
  theme_bw()

# Assign to paths_sf
# st_crs(path_line_sf) <- reef_crs
# st_crs(path_smooth) <- reef_crs

path_line_sf <- st_transform(path_line_sf, 4326)
path_smooth <- st_transform(path_smooth, 4326)
path_smooth.next <- st_transform(path_smooth.next, 4326)

pts_combined <- do.call(rbind, pts_snapped)
pts_sf <- st_as_sf(pts_combined)
pts_sf <- st_transform(pts_sf, 4326)
pts_sf.next <- st_transform(pts_sf.next, 4326)

m <- leaflet() |>
  addProviderTiles(providers$Esri.WorldImagery, group = "Esri Satellite") |>
  addProviderTiles(providers$Esri.OceanBasemap, group = "Esri Ocean") |>
  # addProviderTiles(providers$CartoDB.Positron,  group = "Light") |>
  addRasterImage(
    coralSDM_crop,
    colors  = pal_coral,
    opacity = 0.7,
    project = TRUE,
    group   = "CoralSDM"
  ) |>
  addRasterImage(
    cond_r,
    opacity = 0.7,
    project = TRUE,
    group   = "Conductance"
  ) |>
    addPolylines(
    data  = path_line_sf,
    color = "orange",
    weight = 3,
    group = "Tow Path"
  ) |>
  addPolylines(
    data  = path_smooth,
    color = "springgreen",
    weight = 3,
    group = "Tow Path - Smooth"
  ) |>
  addPolylines(
    data  = path_smooth.next,
    color = "skyblue",
    weight = 3,
    group = "Tow Path - Smooth 2nd"
  ) |>
  addCircleMarkers(
    data = pts_sf,
    radius = 12,
    stroke = TRUE,
    color = "black",
    fill = TRUE,
    fillColor = "white",
    fillOpacity = 1,
    label = ~rownames(pts_sf),
    labelOptions = labelOptions(
      noHide = TRUE, direction = "center",textsize = 25,
      textOnly = TRUE),
    group = "Snapped Points"
  ) |>
  addCircleMarkers(
    data = pts_sf.next,
    radius = 12,
    stroke = TRUE,
    color = "blue",
    fill = TRUE,
    fillColor = "white",
    fillOpacity = 1,
    label = ~rownames(pts_sf.next),
    labelOptions = labelOptions(
      noHide = TRUE, direction = "center",textsize = 25,
      textOnly = TRUE),
    group = "Snapped Points 2nd"
  ) |>
  addLegend(
    pal = pal_coral, values = rng, title = "Probability",
    position = "bottomright", opacity = 1
  ) |>
  addLayersControl(
  baseGroups    = c("Esri Satellite", "Esri Ocean"),
  overlayGroups = c("CoralSDM", "Conductance", "Tow Path", 
                    "Tow Path - Smooth", "Snapped Points", 
                    "Snapped Points 2nd", "Tow Path - Smooth 2nd"),
  options = layersControlOptions(collapsed = FALSE)
)|>
  leafem::addMouseCoordinates()

m